########################################################################################################################
##producing sub_Patristic distance matrixes for pairs of islands

myPatrisDis=read.table("CO1PatrisDist.txt",header=TRUE)
thres=load("list of island3to4.RData")

#########
aux_pdist=NULL
aux_matrix=NULL
PatrDis_y=NULL
for (s in 1:length(thres))
	{
	aux_matrix=thres[[s]]
	aux_pdist=row.names(myPatrisDis) %in% aux_matrix[,1]
	PatrDis_y=subset(myPatrisDis,aux_pdist,aux_pdist)
	
	###sub_matrix with only Moorea (m) and Tahiti (t) islands Patristic distances
	aux_PatrDist_mt=row.names(PatrDis_y) %in% c(grep("moor",row.names(PatrDis_y),fixed=TRUE,value=TRUE),grep("tahi",row.names(PatrDis_y),fixed=TRUE,value=TRUE))
	PatrDist_mt=subset(PatrDis_y,aux_PatrDist_mt,aux_PatrDist_mt)
	write.table(PatrDist_mt,paste("patristic distances Moorea_Tahiti_",names(thres[s]),".txt"))
	
	###sub_matrix with only Moorea (m) and Huahine (h) islands Patristic distances
	aux_PatrDist_mh=row.names(PatrDis_y) %in% c(grep("moor",row.names(PatrDis_y),fixed=TRUE,value=TRUE),grep("huah",row.names(PatrDis_y),fixed=TRUE,value=TRUE))
	PatrDist_mh=subset(PatrDis_y,aux_PatrDist_mh,aux_PatrDist_mh)
	write.table(PatrDist_mh,paste("patristic distances Moorea_Huahine_",names(thres[s]),".txt"))
	
	###sub_matrix with only Moorea (m) and Raiatea (h) islands Patristic distances
	aux_PatrDist_mr=row.names(PatrDis_y) %in% c(grep("moor",row.names(PatrDis_y),fixed=TRUE,value=TRUE),grep("raia",row.names(PatrDis_y),fixed=TRUE,value=TRUE))
	PatrDist_mr=subset(PatrDis_y,aux_PatrDist_mr,aux_PatrDist_mr)
	write.table(PatrDist_mr,paste("patristic distances Moorea_Raiatea_",names(thres[s]),".txt"))
	
	###sub_matrix with only Huahine (m) and Raiatea (h) islands Patristic distances
	aux_PatrDist_hr=row.names(PatrDis_y) %in% c(grep("huah",row.names(PatrDis_y),fixed=TRUE,value=TRUE),grep("raia",row.names(PatrDis_y),fixed=TRUE,value=TRUE))
	PatrDist_hr=subset(PatrDis_y,aux_PatrDist_hr,aux_PatrDist_hr)
	write.table(PatrDist_hr,paste("patristic distances Huahine_Raiatea_",names(thres[s]),".txt"))
	
	###sub_matrix with only Raiatea (m) and Tahiti (h) islands Patristic distances
	aux_PatrDist_rt=row.names(PatrDis_y) %in% c(grep("raia",row.names(PatrDis_y),fixed=TRUE,value=TRUE),grep("tahi",row.names(PatrDis_y),fixed=TRUE,value=TRUE))
	PatrDist_rt=subset(PatrDis_y,aux_PatrDist_rt,aux_PatrDist_rt)
	write.table(PatrDist_rt,paste("patristic distances Raiatea_Tahiti_",names(thres[s]),".txt"))
	
	###sub_matrix with only Huahine (m) and Tahiti (h) islands Patristic distances
	aux_PatrDist_ht=row.names(PatrDis_y) %in% c(grep("huah",row.names(PatrDis_y),fixed=TRUE,value=TRUE),grep("tahi",row.names(PatrDis_y),fixed=TRUE,value=TRUE))
	PatrDist_ht=subset(PatrDis_y,aux_PatrDist_ht,aux_PatrDist_ht)
	write.table(PatrDist_ht,paste("patristic distances Huahine_Tahiti_",names(thres[s]),".txt"))


##producing histograms with frequency of CO1 patristic distances between pairs of islands

mt_vector=as.vector(as.matrix(PatrDist_mt),mode="numeric")
mr_vector=as.vector(as.matrix(PatrDist_mr),mode="numeric")
mh_vector=as.vector(as.matrix(PatrDist_mh),mode="numeric")
hr_vector=as.vector(as.matrix(PatrDist_hr),mode="numeric")
ht_vector=as.vector(as.matrix(PatrDist_ht),mode="numeric")
rt_vector=as.vector(as.matrix(PatrDist_rt),mode="numeric")

pdf(paste("Histogram CO1 Distance Moorea-Tahiti ",names(thres[s]),".pdf"))
hist_mt=hist(mt_vector, breaks=130,col="deepskyblue",xlim=c(0,3), border="royalblue2",main=paste("Histogram for Moorea-Tahiti",names(thres[s])),xlab="CO1 patristic distances",ylab="counts",cex=0.8)
dev.off()

pdf(paste("Histogram CO1 Distance Moorea-Raiatea",names(thres[s]),".pdf"))
hist_mr=hist(mr_vector, breaks=130,col="deepskyblue",xlim=c(0,3), border="royalblue2",main=paste("Histogram for Moorea-Raiatea", names(thres[s])),xlab="CO1 patristic distances",ylab="counts",cex=0.8)
dev.off()

pdf(paste("Histogram CO1 Distance Moorea-Huahine",names(thres[s]),".pdf"))
hist_mh=hist(mh_vector, breaks=125,col="deepskyblue",xlim=c(0,3), border="royalblue2",main=paste("Histogram for Moorea-Huahine",names(thres[s])),xlab="CO1 patristic distances",ylab="counts",cex=0.8)
dev.off()

pdf(paste("Histogram CO1 Distance Huahine-Raiatea",names(thres[s]),".pdf"))
hist_hr=hist(hr_vector, breaks=125,col="deepskyblue",xlim=c(0,3), border="royalblue2",main=paste("Histogram for Huahine-Raiatea",names(thres[s])),xlab="CO1 patristic distances",ylab="counts",cex=0.8)
dev.off()

pdf(paste("Histogram CO1 Distance Huahine-Tahiti",names(thres[s]),".pdf"))
hist_ht=hist(ht_vector, breaks=125,col="deepskyblue",xlim=c(0,3), border="royalblue2",main=paste("Histogram for Huahine-Tahiti",names(thres[s])),xlab="CO1 patristic distances",ylab="counts",cex=0.8)
dev.off()

pdf(paste("Histogram CO1 Distance Raiatea-Tahiti",names(thres[s]),".pdf"))
hist_rt=hist(rt_vector, breaks=125,col="deepskyblue",xlim=c(0,3), border="royalblue2",main=paste("Histogram for Raiatea-Tahiti",names(thres[s])),xlab="CO1 patristic distances",ylab="counts",cex=0.8)
dev.off()
}